R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## report mean ki per ROI

dataMNI <- read.csv('~/Dropbox/DCCN/Creativity/roi_results_mni.csv')
  dataMNI %>%
  group_by(ROI) %>%
  dplyr::summarise(mean = mean(meanRoiContrastEstimate),
  sd = sd(meanRoiContrastEstimate)) 
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ stringr 1.4.0
## ✓ tidyr   1.1.2     ✓ forcats 0.5.0
## ✓ readr   1.4.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## graph the distribution in box plot

ggplot(dataMNI, aes(x=ROI, y=meanRoiContrastEstimate)) + 
    geom_boxplot() + 
    coord_flip()

## look at the density plot

ggplot(dataMNI, aes(x=meanRoiContrastEstimate, group=ROI)) + 
    geom_line(stat="density")

## report mean ki per ROI

dataNative <- read.csv('~/Dropbox/DCCN/Creativity/roi_results_native.csv')
  dataNative %>%
  group_by(ROI) %>%
  dplyr::summarise(mean = mean(meanRoiContrastEstimate),
  sd = sd(meanRoiContrastEstimate)) 
## graph the distribution in box plot

ggplot(dataNative, aes(x=ROI, y=meanRoiContrastEstimate)) + 
    geom_boxplot() + 
    coord_flip()

## look at the density plot

ggplot(dataNative, aes(x=meanRoiContrastEstimate, group=ROI)) + 
    geom_line(stat="density")

## Ruben shared that we are only using native space and only whole caudate, putamen and VS
library(sjPlot)
dataAll <- read.csv('~/Dropbox/DCCN/Creativity/CreativityIDsheet_MissingDataRemoved.csv')

# turn drug conditions into factor levels
dataAll$Drug <- factor(dataAll$Drug, levels = c("MPH","SUL","PBO"));

# set contrasts to sum-to-zero
options(contrasts=c("contr.sum", "contr.poly"))

# set two dataframes for the contrast between MPH and PBO - between SUL and PBO
df_MPH <- dplyr::filter(dataAll, Drug %in% c("MPH","PBO"))  
df_SUL <- dplyr::filter(dataAll, Drug %in% c("SUL","PBO"))
SelectData <- dataNative[grepl("wholeCaudate", dataNative[["ROI"]]) | grepl("wholePutamen", dataNative[["ROI"]]) | grepl("VS", dataNative[["ROI"]]), ]

SelectData <- SelectData[-grep("_", SelectData$ROI),]

ggplot(SelectData, aes(x=ROI, y=meanRoiContrastEstimate)) + 
    geom_boxplot() + 
    coord_flip()

ggplot(SelectData, aes(x=meanRoiContrastEstimate, group=ROI)) + 
    geom_line(stat="density")

# How do different ROI Ki s relate to each other?
my_cols <- c("#00AFBB", "#E7B800", "#FC4E07")  
pairs(dataAll[,2:4], pch = 19,  cex = 0.5,
      lower.panel=NULL)

library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
# Does session number have an effect on convergent thinking?

ggplot(dataAll, aes(x=as.factor(Session), y=Convergent_Pasta)) + 
    geom_boxplot(fill="slateblue", alpha=0.2) + 
    xlab("Session")

# Total number of generated ideas

SessionEffectTotal <- lmer(Total ~ 1 + Session + (1 + Session|| ID), data = dataAll, REML=F) 
print(summary(SessionEffectTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Total ~ 1 + Session + ((1 | ID) + (0 + Session | ID))
##    Data: dataAll
## 
##      AIC      BIC   logLik deviance df.resid 
##   2202.0   2220.2  -1096.0   2192.0      277 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0605 -0.4235 -0.0551  0.3522  4.1868 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 167.328  12.936  
##  ID.1     Session       9.404   3.067  
##  Residual              56.000   7.483  
## Number of obs: 282, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  28.0887     1.7805  15.776
## Session       2.3457     0.6308   3.719
plot_model(SessionEffectTotal, type = "pred", show.data=TRUE, terms = c("Session"))

# Pure convergent scores

SessionEffect <- lmer(Convergent_Pasta ~ 1 + Session + (1 + Session|| ID), data = dataAll, REML=F) 
print(summary(SessionEffect), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Convergent_Pasta ~ 1 + Session + ((1 | ID) + (0 + Session | ID))
##    Data: dataAll
## 
##      AIC      BIC   logLik deviance df.resid 
##   2183.2   2201.4  -1086.6   2173.2      277 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1205 -0.4193 -0.0678  0.3707  5.0262 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 144.02   12.001  
##  ID.1     Session       9.27    3.045  
##  Residual              53.67    7.326  
## Number of obs: 282, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  16.9078     1.6924   9.990
## Session       2.8777     0.6198   4.643
plot_model(SessionEffect, type = "pred", show.data=TRUE, terms = c("Session"))

# convergent/divergent ratio

SessionEffectRatio <- lmer(Con_Div ~ 1 + Session + (1 + Session|| ID), data = dataAll, REML=F) 
print(summary(SessionEffectRatio), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Div ~ 1 + Session + ((1 | ID) + (0 + Session | ID))
##    Data: dataAll
## 
##      AIC      BIC   logLik deviance df.resid 
##   1796.2   1814.1   -893.1   1786.2      260 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0127 -0.2844 -0.2047 -0.0056  5.2673 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept)  9.917   3.149   
##  ID.1     Session      3.854   1.963   
##  Residual             31.422   5.606   
## Number of obs: 265, groups:  ID, 93
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   3.0677     0.9668   3.173
## Session       0.9149     0.4753   1.925
plot_model(SessionEffectRatio, type = "pred", show.data=TRUE,terms = c("Session"))

# Does session number have an effect on divergent thinking?

ggplot(dataAll, aes(x=as.factor(Session), y=Divergent_Pasta)) + 
    geom_boxplot(fill="slateblue", alpha=0.2) + 
    xlab("Session")

SessionEffect2 <- lmer(Divergent_Pasta ~ 1 + Session + (1 + Session|| ID), data = dataAll, REML=F) 
print(summary(SessionEffect2), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Divergent_Pasta ~ 1 + Session + ((1 | ID) + (0 + Session | ID))
##    Data: dataAll
## 
##      AIC      BIC   logLik deviance df.resid 
##   1872.3   1890.5   -931.1   1862.3      277 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2533 -0.5594 -0.1156  0.4431  3.5352 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 27.342   5.229   
##  ID.1     Session      2.999   1.732   
##  Residual             21.882   4.678   
## Number of obs: 282, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  11.1809     0.9133  12.243
## Session      -0.5319     0.3851  -1.381
plot_model(SessionEffect2, type = "pred", show.data=TRUE, terms = c("Session"))

## look at the correlation between convergent and divergent thinking for each session
my_cols <- c("#00AFBB", "#E7B800", "#FC4E07")  
pairs(dataAll[,7:8], pch = 19,  cex = 0.5,
      col = my_cols[dataAll$Session],
      lower.panel=NULL)

# Does drug on its own have an effect on Convergent thinking?

ggplot(dataAll, aes(x=as.factor(Drug), y=Convergent_Pasta)) + 
    geom_boxplot(fill="slateblue", alpha=0.2) + 
    xlab("Session")

DrugEffect <- lmer(Convergent_Pasta ~ 1 + Drug + (1 | ID), data = df_MPH, REML=F) 
print(summary(DrugEffect), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Convergent_Pasta ~ 1 + Drug + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1523.8   1536.8   -757.9   1515.8      184 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8151 -0.4594 -0.0832  0.3814  3.7880 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 159.59   12.633  
##  Residual              85.41    9.242  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  22.7447     1.4670  15.504
## Drug1         0.3085     0.6740   0.458
plot_model(DrugEffect, type = "pred", show.data=TRUE, terms = c("Drug"))

# Drug effect in the presence of session?

DrugAndSessionEffect <- lmer(Convergent_Pasta ~ 1 + Drug*Session + (1 | ID), data = df_MPH, REML=F) 
print(summary(DrugAndSessionEffect), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Convergent_Pasta ~ 1 + Drug * Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1520.3   1539.7   -754.1   1508.3      182 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2378 -0.4158 -0.0723  0.3572  4.1477 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 165.74   12.874  
##  Residual              77.86    8.824  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    18.1707     2.2385   8.117
## Drug1          -1.8587     2.5543  -0.728
## Session         2.3334     0.8697   2.683
## Drug1:Session   1.0552     1.2772   0.826
plot_model(DrugAndSessionEffect, type = "pred", show.data=TRUE, terms = c("Drug","Session"))

# Does drug on its own have an effect on Divergent thinking?

ggplot(dataAll, aes(x=as.factor(Drug), y=Divergent_Pasta)) + 
    geom_boxplot(fill="slateblue", alpha=0.2) + 
    xlab("Session")

DrugEffect2 <- lmer(Divergent_Pasta ~ 1 + Drug + (1 | ID), data = df_MPH, REML=F) 
print(summary(DrugEffect2), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Divergent_Pasta ~ 1 + Drug + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1290.9   1303.8   -641.4   1282.9      184 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3879 -0.4947 -0.0901  0.3593  3.6678 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 37.76    6.145   
##  Residual             28.00    5.292   
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  10.1862     0.7421  13.727
## Drug1         0.6968     0.3859   1.805
plot_model(DrugEffect2, type = "pred", terms = c("Drug"))

# Drug effect in the presence of session?

DrugAndSessionEffect2 <- lmer(Divergent_Pasta ~ 1 + Drug* Session + (1 | ID), data = df_MPH, REML=F) 
print(summary(DrugAndSessionEffect2), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Divergent_Pasta ~ 1 + Drug * Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1292.0   1311.4   -640.0   1280.0      182 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2919 -0.4825 -0.0912  0.3654  3.7287 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 36.27    6.023   
##  Residual             27.97    5.289   
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    11.1113     1.2360   8.990
## Drug1           2.6937     1.4175   1.900
## Session        -0.4498     0.5148  -0.874
## Drug1:Session  -1.0190     0.7047  -1.446
plot_model(DrugAndSessionEffect2, type = "pred", terms = c("Drug","Session"))

# How do different ROI Ki s relate to each other?
cols <- c(2:4, 7:8)
my_cols <- c("#00AFBB", "#E7B800", "#FC4E07")  
pairs(dataAll[,cols], pch = 19,  cex = 0.5,
            col = my_cols[dataAll$Session],
      lower.panel=NULL)

############################# Test the effects of Caudate Ki #############################################################

## effects of MPH vs PBO below

# pure Convergent score
MPH_Caudate_PureConvergent <- lmer(Convergent_Pasta ~ 1 + Drug*Caudate_ki + Session + (1 | ID), data = df_MPH, REML=F) 
print(summary(MPH_Caudate_PureConvergent), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Convergent_Pasta ~ 1 + Drug * Caudate_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1522.0   1544.7   -754.0   1508.0      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1830 -0.4290 -0.0762  0.3575  4.1319 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 165.85   12.878  
##  Residual              77.65    8.812  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       18.7640    10.3977   1.805
## Drug1             -4.0905     4.4653  -0.916
## Caudate_ki       -50.2712   722.7969  -0.070
## Session            2.4191     0.8713   2.776
## Drug1:Caudate_ki 305.3485   315.6846   0.967
plot_model(MPH_Caudate_PureConvergent, type = "pred", terms = c("Caudate_ki","Session","Drug"))

plot(MPH_Caudate_PureConvergent)  # Plot the model information

qqnorm(residuals(MPH_Caudate_PureConvergent)) 

# Convergent/Divergent ratio score
MPH_Caudate_ConDivRatio <- lmer(Con_Div ~ 1 + Drug*Caudate_ki + Session + (1 | ID), data = df_MPH, REML=F)
print(summary(MPH_Caudate_ConDivRatio), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Div ~ 1 + Drug * Caudate_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1226.7   1248.8   -606.3   1212.7      168 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7803 -0.3750 -0.2319 -0.0082  5.1290 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 20.16    4.490   
##  Residual             42.87    6.547   
## Number of obs: 175, groups:  ID, 91
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)        3.6777     4.8746   0.754
## Drug1             -3.5583     3.3911  -1.049
## Caudate_ki       -27.9163   330.1346  -0.085
## Session            0.8616     0.6456   1.335
## Drug1:Caudate_ki 253.1895   239.3929   1.058
plot_model(MPH_Caudate_ConDivRatio, type = "pred", terms = c("Caudate_ki","Session","Drug"))

plot(MPH_Caudate_ConDivRatio)

qqnorm(residuals(MPH_Caudate_ConDivRatio))

# Convergent/Total ratio score
MPH_Caudate_ConTotal <- lmer(Con_Total ~ 1 + Drug*Caudate_ki + Session + (1 | ID), data = df_MPH, REML=F)
print(summary(MPH_Caudate_ConTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Total ~ 1 + Drug * Caudate_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##    -54.5    -31.9     34.3    -68.5      179 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0574 -0.5147  0.1010  0.5240  2.1860 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.03478  0.1865  
##  Residual             0.01852  0.1361  
## Number of obs: 186, groups:  ID, 94
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       0.50870    0.15288   3.327
## Drug1            -0.05650    0.06897  -0.819
## Caudate_ki        6.22253   10.60444   0.587
## Session           0.03407    0.01357   2.511
## Drug1:Caudate_ki  3.18761    4.87762   0.654
plot_model(MPH_Caudate_ConTotal, type = "pred", terms = c("Caudate_ki","Session","Drug"))

plot(MPH_Caudate_ConTotal)

qqnorm(residuals(MPH_Caudate_ConTotal))

# Convergent-Divergent difference score 
MPH_Caudate_ConDifference <- lmer(ConDiv_Dif ~ 1 + Drug*Caudate_ki + Session + (1 | ID), data = df_MPH, REML=F)
print(summary(MPH_Caudate_ConDifference), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: ConDiv_Dif ~ 1 + Drug * Caudate_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1609.7   1632.4   -797.9   1595.7      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1195 -0.4122 -0.0729  0.4347  3.9454 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 232.5    15.25   
##  Residual             134.8    11.61   
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)         6.365     12.624   0.504
## Drug1              -3.333      5.883  -0.566
## Caudate_ki         56.472    874.925   0.065
## Session             2.791      1.142   2.443
## Drug1:Caudate_ki  199.909    415.911   0.481
plot_model(MPH_Caudate_ConDifference, type = "pred", terms = c("Caudate_ki","Session","Drug"))

plot(MPH_Caudate_ConDifference)

qqnorm(residuals(MPH_Caudate_ConDifference))

# Divergent/Total ratio score
MPH_Caudate_DivTotal <- lmer(Div_Total ~ 1 + Drug*Caudate_ki + Session + (1 | ID), data = df_MPH, REML=F)
print(summary(MPH_Caudate_DivTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Div_Total ~ 1 + Drug * Caudate_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##    -54.4    -31.8     34.2    -68.4      179 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1847 -0.5250 -0.1048  0.5143  2.0560 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.03477  0.1865  
##  Residual             0.01855  0.1362  
## Number of obs: 186, groups:  ID, 94
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       0.49187    0.15288   3.217
## Drug1             0.05697    0.06903   0.825
## Caudate_ki       -6.22980   10.60419  -0.587
## Session          -0.03418    0.01358  -2.517
## Drug1:Caudate_ki -3.21665    4.88194  -0.659
plot_model(MPH_Caudate_DivTotal, type = "pred", terms = c("Caudate_ki","Session","Drug"))

plot(MPH_Caudate_DivTotal)

qqnorm(residuals(MPH_Caudate_DivTotal))

## effects of SUL vs PBO below

# pure Convergent score
SUL_Caudate_PureConvergent <- lmer(Convergent_Pasta ~ 1 + Drug*Caudate_ki + Session + (1 | ID), data = df_SUL, REML=F) 
print(summary(SUL_Caudate_PureConvergent), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Convergent_Pasta ~ 1 + Drug * Caudate_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##   1494.8   1517.5   -740.4   1480.8      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9468 -0.4408 -0.0731  0.2915  3.9815 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 178.63   13.365  
##  Residual              57.41    7.577  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       18.1860    10.3607   1.755
## Drug1             -4.0681     3.8259  -1.063
## Caudate_ki       -88.7444   727.4964  -0.122
## Session            2.7543     0.7906   3.484
## Drug1:Caudate_ki 269.0754   270.9883   0.993
plot_model(SUL_Caudate_PureConvergent, type = "pred", terms = c("Caudate_ki","Session","Drug"))

plot(SUL_Caudate_PureConvergent)

qqnorm(residuals(SUL_Caudate_PureConvergent))

# Convergent/Divergent ratio score
SUL_Caudate_ConDivRatio <- lmer(Con_Div ~ 1 + Drug*Caudate_ki + Session + (1 | ID), data = df_SUL, REML=F)
print(summary(SUL_Caudate_ConDivRatio), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Div ~ 1 + Drug * Caudate_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##   1206.4   1228.7   -596.2   1192.4      171 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8341 -0.3568 -0.2152  0.0065  5.7767 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 16.10    4.012   
##  Residual             33.96    5.827   
## Number of obs: 178, groups:  ID, 93
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)         8.3899     4.2775   1.961
## Drug1               0.4849     3.0083   0.161
## Caudate_ki       -369.9599   294.5373  -1.256
## Session             0.6971     0.5930   1.175
## Drug1:Caudate_ki  -61.0075   213.5831  -0.286
plot_model(SUL_Caudate_ConDivRatio, type = "pred", terms = c("Caudate_ki","Session","Drug"))

plot(SUL_Caudate_ConDivRatio)

qqnorm(residuals(SUL_Caudate_ConDivRatio))

# Convergent/Total ratio score
SUL_Caudate_ConTotal <- lmer(Con_Total ~ 1 + Drug*Caudate_ki + Session + (1 | ID), data = df_SUL, REML=F)
print(summary(SUL_Caudate_ConTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Total ~ 1 + Drug * Caudate_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##    -50.9    -28.2     32.4    -64.9      181 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.22401 -0.55586  0.07603  0.51605  1.83188 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.03359  0.1833  
##  Residual             0.01977  0.1406  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)       0.497809   0.150729   3.303
## Drug1            -0.006216   0.071004  -0.088
## Caudate_ki        1.950943  10.536689   0.185
## Session           0.068229   0.014451   4.721
## Drug1:Caudate_ki -0.859790   5.029105  -0.171
plot_model(SUL_Caudate_ConTotal, type = "pred", terms = c("Caudate_ki","Session","Drug"))

plot(SUL_Caudate_ConTotal)

qqnorm(residuals(SUL_Caudate_ConTotal))

# Convergent-Divergent difference score 
SUL_Caudate_ConDifference <- lmer(ConDiv_Dif ~ 1 + Drug*Caudate_ki + Session + (1 | ID), data = df_SUL, REML=F)
print(summary(SUL_Caudate_ConDifference), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: ConDiv_Dif ~ 1 + Drug * Caudate_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##   1589.6   1612.3   -787.8   1575.6      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7010 -0.4469 -0.0126  0.3783  3.7805 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 217.3    14.74   
##  Residual             118.1    10.87   
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)         5.339     12.009   0.445
## Drug1              -3.695      5.487  -0.673
## Caudate_ki         75.444    840.056   0.090
## Session             3.162      1.119   2.825
## Drug1:Caudate_ki  221.314    388.615   0.569
plot_model(SUL_Caudate_ConDifference, type = "pred", terms = c("Caudate_ki","Session","Drug"))

plot(SUL_Caudate_ConDifference)

qqnorm(residuals(SUL_Caudate_ConDifference))

# Divergent/Total ratio score
SUL_Caudate_DivTotal <- lmer(Div_Total ~ 1 + Drug*Caudate_ki + Session + (1 | ID), data = df_SUL, REML=F)
print(summary(SUL_Caudate_DivTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Div_Total ~ 1 + Drug * Caudate_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##    -50.9    -28.3     32.5    -64.9      181 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.83278 -0.51699 -0.07626  0.55621  2.22345 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.03356  0.1832  
##  Residual             0.01977  0.1406  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)       0.502142   0.150691   3.332
## Drug1             0.006258   0.071005   0.088
## Caudate_ki       -1.940537  10.533992  -0.184
## Session          -0.068224   0.014451  -4.721
## Drug1:Caudate_ki  0.849165   5.029166   0.169
plot_model(SUL_Caudate_DivTotal, type = "pred", terms = c("Caudate_ki","Session","Drug"))

plot(SUL_Caudate_DivTotal)

qqnorm(residuals(SUL_Caudate_DivTotal))

############################# Test the effects of Putamen Ki #############################################################

## effects of MPH vs PBO below

# pure Convergent score
MPH_Putamen_PureConvergent <- lmer(Convergent_Pasta ~ 1 + Drug*Putamen_ki + Session + (1 | ID), data = df_MPH, REML=F) 
print(summary(MPH_Putamen_PureConvergent), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Convergent_Pasta ~ 1 + Drug * Putamen_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1522.9   1545.6   -754.5   1508.9      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1949 -0.4529 -0.0786  0.3663  4.1514 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 165.35   12.859  
##  Residual              78.46    8.858  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)       19.36174   11.43444   1.693
## Drug1              0.04279    5.02649   0.009
## Putamen_ki       -69.91372  664.58081  -0.105
## Session            2.35409    0.88805   2.651
## Drug1:Putamen_ki   8.36032  296.07382   0.028
plot_model(MPH_Putamen_PureConvergent, type = "pred", terms = c("Putamen_ki","Session","Drug"))

# Convergent/Divergent ratio score
MPH_Putamen_ConDivRatio <- lmer(Con_Div ~ 1 + Drug*Putamen_ki + Session + (1 | ID), data = df_MPH, REML=F)
print(summary(MPH_Putamen_ConDivRatio), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Div ~ 1 + Drug * Putamen_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1227.6   1249.8   -606.8   1213.6      168 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8090 -0.3680 -0.2338 -0.0256  5.1362 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 20.20    4.495   
##  Residual             43.15    6.569   
## Number of obs: 175, groups:  ID, 91
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)        4.1489     5.3267   0.779
## Drug1              1.4390     3.7759   0.381
## Putamen_ki       -40.7150   303.5188  -0.134
## Session            0.7783     0.6553   1.188
## Drug1:Putamen_ki -86.0849   222.2221  -0.387
plot_model(MPH_Putamen_ConDivRatio, type = "pred", terms = c("Putamen_ki","Session","Drug"))

# Convergent/Total ratio score
MPH_Putamen_ConTotal <- lmer(Con_Total ~ 1 + Drug*Putamen_ki + Session + (1 | ID), data = df_MPH, REML=F)
print(summary(MPH_Putamen_ConTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Total ~ 1 + Drug * Putamen_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##    -54.1    -31.5     34.0    -68.1      179 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0773 -0.4750  0.0937  0.5361  2.1925 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.03480  0.1865  
##  Residual             0.01858  0.1363  
## Number of obs: 186, groups:  ID, 94
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)       0.510966   0.168119   3.039
## Drug1             0.006385   0.077366   0.083
## Putamen_ki        5.199394   9.757426   0.533
## Session           0.032703   0.013777   2.374
## Drug1:Putamen_ki -1.086851   4.556358  -0.239
plot_model(MPH_Putamen_ConTotal, type = "pred", terms = c("Putamen_ki","Session","Drug"))

# Convergent-Divergent difference score 
MPH_Putamen_ConDifference <- lmer(ConDiv_Dif ~ 1 + Drug*Putamen_ki + Session + (1 | ID), data = df_MPH, REML=F)
print(summary(MPH_Putamen_ConDifference), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: ConDiv_Dif ~ 1 + Drug * Putamen_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1609.9   1632.5   -797.9   1595.9      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1311 -0.3807 -0.0473  0.4687  3.9547 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 232.2    15.24   
##  Residual             135.0    11.62   
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)         9.432     13.874   0.680
## Drug1               1.538      6.593   0.233
## Putamen_ki       -121.965    804.469  -0.152
## Session             2.674      1.159   2.306
## Drug1:Putamen_ki -123.085    388.327  -0.317
plot_model(MPH_Putamen_ConDifference, type = "pred", terms = c("Putamen_ki","Session","Drug"))

plot(MPH_Putamen_ConDifference)

qqnorm(residuals(MPH_Putamen_ConDifference))

# Divergent/Total ratio score
MPH_Putamen_DivTotal <- lmer(Div_Total ~ 1 + Drug*Putamen_ki + Session + (1 | ID), data = df_MPH, REML=F)
print(summary(MPH_Putamen_DivTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Div_Total ~ 1 + Drug * Putamen_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##    -53.9    -31.3     34.0    -67.9      179 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.19125 -0.53605 -0.09337  0.47557  2.07518 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.03478  0.1865  
##  Residual             0.01862  0.1364  
## Number of obs: 186, groups:  ID, 94
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)       0.489172   0.168122   2.910
## Drug1            -0.006116   0.077438  -0.079
## Putamen_ki       -5.179289   9.757352  -0.531
## Session          -0.032811   0.013789  -2.380
## Drug1:Putamen_ki  1.074502   4.560597   0.236
plot_model(MPH_Putamen_DivTotal, type = "pred", terms = c("Putamen_ki","Session","Drug"))

## effects of SUL vs PBO below

# pure Convergent score
SUL_Putamen_PureConvergent <- lmer(Convergent_Pasta ~ 1 + Drug*Putamen_ki + Session + (1 | ID), data = df_SUL, REML=F) 
print(summary(SUL_Putamen_PureConvergent), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Convergent_Pasta ~ 1 + Drug * Putamen_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##   1495.8   1518.4   -740.9   1481.8      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9680 -0.4589 -0.1037  0.3027  4.0126 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 178.33   13.354  
##  Residual              57.99    7.615  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       17.6047    11.3887   1.546
## Drug1             -1.2096     4.2438  -0.285
## Putamen_ki       -44.8624   669.5883  -0.067
## Session            2.8012     0.7939   3.528
## Drug1:Putamen_ki  53.4798   250.1742   0.214
plot_model(SUL_Putamen_PureConvergent, type = "pred", terms = c("Putamen_ki","Session","Drug"))

# Convergent/Divergent ratio score
SUL_Putamen_ConDivRatio <- lmer(Con_Div ~ 1 + Drug*Putamen_ki + Session + (1 | ID), data = df_SUL, REML=F)
print(summary(SUL_Putamen_ConDivRatio), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Div ~ 1 + Drug * Putamen_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##   1205.6   1227.9   -595.8   1191.6      171 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8730 -0.3730 -0.2144  0.0592  5.7171 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 17.04    4.128   
##  Residual             33.11    5.754   
## Number of obs: 178, groups:  ID, 93
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)         7.1114     4.6993   1.513
## Drug1               4.0429     3.3049   1.223
## Putamen_ki       -232.5022   274.6479  -0.847
## Session             0.7103     0.5885   1.207
## Drug1:Putamen_ki -262.3971   194.9995  -1.346
plot_model(SUL_Putamen_ConDivRatio, type = "pred", terms = c("Putamen_ki","Session","Drug"))

# Convergent/Total ratio score
SUL_Putamen_ConTotal <- lmer(Con_Total ~ 1 + Drug*Putamen_ki + Session + (1 | ID), data = df_SUL, REML=F)
print(summary(SUL_Putamen_ConTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Total ~ 1 + Drug * Putamen_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##    -51.1    -28.4     32.5    -65.1      181 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.17868 -0.55058  0.07223  0.50396  1.82385 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.03357  0.1832  
##  Residual             0.01975  0.1405  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       0.47846    0.16528   2.895
## Drug1             0.01354    0.07831   0.173
## Putamen_ki        2.82104    9.69880   0.291
## Session           0.06784    0.01444   4.699
## Drug1:Putamen_ki -1.88813    4.61658  -0.409
plot_model(SUL_Putamen_ConTotal, type = "pred", terms = c("Putamen_ki","Session","Drug"))

# Convergent-Divergent difference score 
SUL_Putamen_ConDifference <- lmer(ConDiv_Dif ~ 1 + Drug*Putamen_ki + Session + (1 | ID), data = df_SUL, REML=F)
print(summary(SUL_Putamen_ConDifference), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: ConDiv_Dif ~ 1 + Drug * Putamen_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##   1589.9   1612.6   -788.0   1575.9      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7219 -0.4398 -0.0209  0.3943  3.8002 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 217.2    14.74   
##  Residual             118.5    10.88   
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)        6.5821    13.1780   0.499
## Drug1             -0.7163     6.0655  -0.118
## Putamen_ki       -15.8827   773.5282  -0.021
## Session            3.2009     1.1208   2.856
## Drug1:Putamen_ki   6.6629   357.5738   0.019
plot_model(SUL_Putamen_ConDifference, type = "pred", terms = c("Putamen_ki","Session","Drug"))

plot(SUL_Putamen_ConDifference)

qqnorm(residuals(SUL_Putamen_ConDifference))

# Divergent/Total ratio score
SUL_Putamen_DivTotal <- lmer(Div_Total ~ 1 + Drug*Putamen_ki + Session + (1 | ID), data = df_SUL, REML=F)
print(summary(SUL_Putamen_DivTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Div_Total ~ 1 + Drug * Putamen_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##    -51.1    -28.4     32.6    -65.1      181 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.82483 -0.50493 -0.07243  0.55023  2.17833 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.03355  0.1832  
##  Residual             0.01975  0.1405  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       0.52141    0.16524   3.156
## Drug1            -0.01341    0.07831  -0.171
## Putamen_ki       -2.80731    9.69635  -0.290
## Session          -0.06784    0.01444  -4.698
## Drug1:Putamen_ki  1.87433    4.61667   0.406
plot_model(SUL_Putamen_DivTotal, type = "pred", terms = c("Putamen_ki","Session","Drug"))

############################# Test the effects of VS Ki #############################################################

## effects of MPH vs PBO below

# pure Convergent score
MPH_VS_PureConvergent <- lmer(Convergent_Pasta ~ 1 + Drug*VS_ki + Session + (1 | ID), data = df_MPH, REML=F) 
print(summary(MPH_VS_PureConvergent), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Convergent_Pasta ~ 1 + Drug * VS_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1520.0   1542.6   -753.0   1506.0      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2164 -0.4128 -0.0839  0.3701  4.1958 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 163.41   12.783  
##  Residual              77.03    8.777  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   4.5639    11.5073   0.397
## Drug1        -6.1717     5.0276  -1.228
## VS_ki       907.3052   772.0248   1.175
## Session       2.5662     0.8773   2.925
## Drug1:VS_ki 435.6364   341.8961   1.274
plot_model(MPH_VS_PureConvergent, show.values=TRUE, value.offset = .3)

plot_model(MPH_VS_PureConvergent,  type = "pred", terms = c("VS_ki","Session","Drug"))

# Convergent/Divergent ratio score
MPH_VS_ConDivRatio <- lmer(Con_Div ~ 1 + Drug*VS_ki + Session + (1 | ID), data = df_MPH, REML=F)
print(summary(MPH_VS_ConDivRatio), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Div ~ 1 + Drug * VS_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1227.6   1249.7   -606.8   1213.6      168 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7999 -0.3856 -0.2262 -0.0557  5.1621 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 20.07    4.480   
##  Residual             43.21    6.574   
## Number of obs: 175, groups:  ID, 91
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   1.7153     5.4730   0.313
## Drug1        -1.5541     3.8669  -0.402
## VS_ki       108.2781   359.4776   0.301
## Session       0.8644     0.6529   1.324
## Drug1:VS_ki 105.5954   262.7893   0.402
plot_model(MPH_VS_ConDivRatio, show.values=TRUE, value.offset = .3)

plot_model(MPH_VS_ConDivRatio, type = "pred", terms = c("VS_ki","Session","Drug"))

# Convergent/Total ratio score
MPH_VS_ConTotal <- lmer(Con_Total ~ 1 + Drug*VS_ki + Session + (1 | ID), data = df_MPH, REML=F)
print(summary(MPH_VS_ConTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Total ~ 1 + Drug * VS_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##    -55.1    -32.5     34.5    -69.1      179 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1040 -0.5145  0.1039  0.5376  2.2060 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.03467  0.1862  
##  Residual             0.01846  0.1359  
## Number of obs: 186, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.47019    0.17019   2.763
## Drug1       -0.08165    0.07822  -1.044
## VS_ki        8.43062   11.39733   0.740
## Session      0.03556    0.01371   2.593
## Drug1:VS_ki  4.79022    5.32597   0.899
plot_model(MPH_VS_ConTotal, show.values=TRUE, value.offset = .3)

plot_model(MPH_VS_ConTotal, type = "pred", terms = c("VS_ki","Session","Drug"))

# Convergent-Divergent difference score 
MPH_VS_ConDifference <- lmer(ConDiv_Dif ~ 1 + Drug*VS_ki + Session + (1 | ID), data = df_MPH, REML=F)
print(summary(MPH_VS_ConDifference), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: ConDiv_Dif ~ 1 + Drug * VS_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##   1608.0   1630.7   -797.0   1594.0      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1317 -0.3960 -0.0366  0.4012  3.9846 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 228.7    15.12   
##  Residual             134.2    11.58   
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   -8.830     13.964  -0.632
## Drug1         -5.837      6.634  -0.880
## VS_ki       1078.450    934.423   1.154
## Session        2.936      1.152   2.549
## Drug1:VS_ki  363.473    451.169   0.806
plot_model(MPH_VS_ConDifference, type = "pred", terms = c("VS_ki","Session","Drug"))

plot(MPH_VS_ConDifference)

qqnorm(residuals(MPH_VS_ConDifference))

# Divergent/Total ratio score
MPH_VS_DivTotal <- lmer(Div_Total ~ 1 + Drug*VS_ki + Session + (1 | ID), data = df_MPH, REML=F)
print(summary(MPH_VS_DivTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Div_Total ~ 1 + Drug * VS_ki + Session + (1 | ID)
##    Data: df_MPH
## 
##      AIC      BIC   logLik deviance df.resid 
##    -54.9    -32.3     34.5    -68.9      179 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2048 -0.5383 -0.1045  0.5148  2.1020 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.03465  0.1861  
##  Residual             0.01850  0.1360  
## Number of obs: 186, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.53043    0.17019   3.117
## Drug1        0.08178    0.07829   1.045
## VS_ki       -8.44068   11.39708  -0.741
## Session     -0.03566    0.01373  -2.598
## Drug1:VS_ki -4.79473    5.33086  -0.899
plot_model(MPH_VS_DivTotal, show.values=TRUE, value.offset = .3)

plot_model(MPH_VS_DivTotal, type = "pred", terms = c("VS_ki","Session","Drug"))

## effects of SUL vs PBO below

# pure Convergent score
SUL_VS_PureConvergent <- lmer(Convergent_Pasta ~ 1 + Drug*VS_ki + Session + (1 | ID), data = df_SUL, REML=F) 
print(summary(SUL_VS_PureConvergent), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Convergent_Pasta ~ 1 + Drug * VS_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##   1493.9   1516.6   -740.0   1479.9      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9575 -0.4389 -0.0915  0.3183  4.0199 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 176.62   13.290  
##  Residual              57.41    7.577  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   5.9742    11.4791   0.520
## Drug1        -4.5189     4.2738  -1.057
## VS_ki       752.8912   779.0124   0.966
## Session       2.7578     0.7898   3.492
## Drug1:VS_ki 289.3428   290.9810   0.994
plot_model(SUL_VS_PureConvergent, show.values=TRUE, value.offset = .3)

plot_model(SUL_VS_PureConvergent, type = "pred", terms = c("VS_ki","Session","Drug"))

# Convergent/Divergent ratio score
SUL_VS_ConDivRatio <- lmer(Con_Div ~ 1 + Drug*VS_ki + Session + (1 | ID), data = df_SUL, REML=F)
print(summary(SUL_VS_ConDivRatio), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Div ~ 1 + Drug * VS_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##   1207.0   1229.3   -596.5   1193.0      171 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8475 -0.3906 -0.1878  0.0142  5.7288 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 16.72    4.089   
##  Residual             33.68    5.803   
## Number of obs: 178, groups:  ID, 93
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)    6.0577     4.8061   1.260
## Drug1          2.4641     3.3765   0.730
## VS_ki       -193.9864   322.5798  -0.601
## Session        0.6976     0.5918   1.179
## Drug1:VS_ki -193.7748   229.5201  -0.844
plot_model(SUL_VS_ConDivRatio, show.values=TRUE, value.offset = .3)

plot_model(SUL_VS_ConDivRatio, type = "pred", terms = c("VS_ki","Session","Drug"))

# Convergent/Total ratio score
SUL_VS_ConTotal <- lmer(Con_Total ~ 1 + Drug*VS_ki + Session + (1 | ID), data = df_SUL, REML=F)
print(summary(SUL_VS_ConTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Con_Total ~ 1 + Drug * VS_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##    -51.3    -28.7     32.7    -65.3      181 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.15283 -0.54150  0.06619  0.52389  1.83639 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.03365  0.1834  
##  Residual             0.01968  0.1403  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.54692    0.16747   3.266
## Drug1        0.03660    0.07912   0.463
## VS_ki       -1.51742   11.34328  -0.134
## Session      0.06837    0.01441   4.744
## Drug1:VS_ki -3.76739    5.38705  -0.699
plot_model(SUL_VS_ConTotal, show.values=TRUE, value.offset = .3)

plot_model(SUL_VS_ConTotal, type = "pred", terms = c("VS_ki","Session","Drug"))

# Convergent-Divergent difference score 
SUL_VS_ConDifference <- lmer(ConDiv_Dif ~ 1 + Drug*VS_ki + Session + (1 | ID), data = df_SUL, REML=F)
print(summary(SUL_VS_ConDifference), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: ConDiv_Dif ~ 1 + Drug * VS_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##   1589.1   1611.8   -787.6   1575.1      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7051 -0.4561 -0.0262  0.3903  3.8065 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 214.9    14.66   
##  Residual             118.4    10.88   
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   -5.174     13.292  -0.389
## Drug1         -1.919      6.138  -0.313
## VS_ki        796.974    900.504   0.885
## Session        3.143      1.120   2.807
## Drug1:VS_ki   90.778    417.924   0.217
plot_model(SUL_VS_ConDifference, type = "pred", terms = c("VS_ki","Session","Drug"))

plot(SUL_VS_ConDifference)

qqnorm(residuals(SUL_VS_ConDifference))

# Divergent/Total ratio score
SUL_VS_DivTotal <- lmer(Div_Total ~ 1 + Drug*VS_ki + Session + (1 | ID), data = df_SUL, REML=F)
print(summary(SUL_VS_DivTotal), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Div_Total ~ 1 + Drug * VS_ki + Session + (1 | ID)
##    Data: df_SUL
## 
##      AIC      BIC   logLik deviance df.resid 
##    -51.4    -28.7     32.7    -65.4      181 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.83744 -0.52377 -0.06621  0.54220  2.15188 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.03362  0.1834  
##  Residual             0.01968  0.1403  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.45326    0.16743   2.707
## Drug1       -0.03678    0.07912  -0.465
## VS_ki        1.51227   11.34036   0.133
## Session     -0.06836    0.01441  -4.744
## Drug1:VS_ki  3.77236    5.38706   0.700
plot_model(SUL_VS_DivTotal, show.values=TRUE, value.offset = .3)

plot_model(SUL_VS_DivTotal, type = "pred", terms = c("VS_ki","Session","Drug"))

# Factorize the data into score types and test the effects of score type as a factor on scores (1=convergent, 2=divergent)

FactorizedScores <- read.csv('~/Dropbox/DCCN/Creativity/CreativityIDsheet_MissingDataRemoved_ScoresFactorized.csv')

# turn drug conditions into factor levels
FactorizedScores$Drug <- factor(FactorizedScores$Drug, levels = c("MPH","SUL","PBO"));
FactorizedScores$ScoreType <- factor(FactorizedScores$ScoreType);
hs <- function(x) {
    y <- 0.5*exp(-x)*(exp(2*x)-1)
    return(y)
}  # inverse sine transormation
#FactorizedScores$Score <- hs(FactorizedScores$Score); # zscored the scores
FactorizedScores$Score <-scale(FactorizedScores$Score, center = TRUE, scale = TRUE)

# ScoreType 1 refers to convergent, 2 refers to divergent

# set contrasts to sum-to-zero
options(contrasts=c("contr.sum", "contr.poly"))

# set two dataframes for the contrast between MPH and PBO - between SUL and PBO
df_MPH_Factorized <- dplyr::filter(FactorizedScores, Drug %in% c("MPH","PBO"))  
df_SUL_Factorized <- dplyr::filter(FactorizedScores, Drug %in% c("SUL","PBO"))
## effects of MPH vs PBO on factorized model

MPH_VS_FactorizedScore <- lmer(Score ~ 1 + Drug*VS_ki*ScoreType + Session + (1 | ID), data = df_MPH_Factorized, REML=F) 
print(summary(MPH_VS_FactorizedScore), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Score ~ 1 + Drug * VS_ki * ScoreType + Session + (1 | ID)
##    Data: df_MPH_Factorized
## 
##      AIC      BIC   logLik deviance df.resid 
##    996.0   1039.2   -487.0    974.0      365 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9628 -0.6380 -0.1269  0.4391  5.9261 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.09562  0.3092  
##  Residual             0.70010  0.8367  
## Number of obs: 376, groups:  ID, 94
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)            -0.46831    0.43371  -1.080
## Drug1                  -0.20533    0.33832  -0.607
## VS_ki                  25.82118   28.28359   0.913
## ScoreType1             -0.07783    0.33367  -0.233
## Session                 0.05047    0.05560   0.908
## Drug1:VS_ki            16.39899   23.01029   0.713
## Drug1:ScoreType1       -0.10384    0.33367  -0.311
## VS_ki:ScoreType1       36.36198   22.72164   1.600
## Drug1:VS_ki:ScoreType1  6.17207   22.72164   0.272
plot_model(MPH_VS_FactorizedScore, show.values=TRUE, value.offset = .3)

plot_model(MPH_VS_FactorizedScore,  type = "pred", terms = c("VS_ki","ScoreType","Drug"))

plot(MPH_VS_FactorizedScore)

qqnorm(residuals(MPH_VS_FactorizedScore))

MPH_Putamen_FactorizedScore <- lmer(Score ~ 1 + Drug*Putamen_ki*ScoreType + Session + (1 | ID), data = df_MPH_Factorized, REML=F) 
print(summary(MPH_Putamen_FactorizedScore), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Score ~ 1 + Drug * Putamen_ki * ScoreType + Session + (1 | ID)
##    Data: df_MPH_Factorized
## 
##      AIC      BIC   logLik deviance df.resid 
##    999.5   1042.7   -488.7    977.5      365 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8945 -0.6345 -0.1318  0.4171  6.0511 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.09574  0.3094  
##  Residual             0.70728  0.8410  
## Number of obs: 376, groups:  ID, 94
## 
## Fixed effects:
##                              Estimate Std. Error t value
## (Intercept)                  -0.05924    0.42951  -0.138
## Drug1                        -0.01969    0.33673  -0.058
## Putamen_ki                   -1.14636   24.26996  -0.047
## ScoreType1                    0.55270    0.33113   1.669
## Session                       0.04335    0.05601   0.774
## Drug1:Putamen_ki              3.18561   19.83745   0.161
## Drug1:ScoreType1              0.16026    0.33113   0.484
## Putamen_ki:ScoreType1        -6.01111   19.53336  -0.308
## Drug1:Putamen_ki:ScoreType1 -10.36694   19.53336  -0.531
plot_model(MPH_Putamen_FactorizedScore, show.values=TRUE, value.offset = .3)

plot_model(MPH_Putamen_FactorizedScore,  type = "pred", terms = c("Putamen_ki","ScoreType","Drug"))

plot(MPH_Putamen_FactorizedScore)

qqnorm(residuals(MPH_Putamen_FactorizedScore))

MPH_Caudate_FactorizedScore <- lmer(Score ~ 1 + Drug*Caudate_ki*ScoreType + Session + (1 | ID), data = df_MPH_Factorized, REML=F) 
print(summary(MPH_Caudate_FactorizedScore), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Score ~ 1 + Drug * Caudate_ki * ScoreType + Session + (1 | ID)
##    Data: df_MPH_Factorized
## 
##      AIC      BIC   logLik deviance df.resid 
##    999.4   1042.6   -488.7    977.4      365 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9110 -0.6321 -0.1090  0.4160  6.0466 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.09567  0.3093  
##  Residual             0.70708  0.8409  
## Number of obs: 376, groups:  ID, 94
## 
## Fixed effects:
##                              Estimate Std. Error t value
## (Intercept)                  0.008465   0.392731   0.022
## Drug1                       -0.160904   0.301161  -0.534
## Caudate_ki                  -6.303045  26.389338  -0.239
## ScoreType1                   0.454306   0.300111   1.514
## Session                      0.043955   0.055316   0.795
## Drug1:Caudate_ki            13.923455  21.293132   0.654
## Drug1:ScoreType1            -0.074264   0.300111  -0.247
## Caudate_ki:ScoreType1       -0.188328  21.233979  -0.009
## Drug1:Caudate_ki:ScoreType1  4.311609  21.233979   0.203
plot_model(MPH_Caudate_FactorizedScore, show.values=TRUE, value.offset = .3)

plot_model(MPH_Caudate_FactorizedScore,  type = "pred", terms = c("Caudate_ki","ScoreType","Drug"))

plot(MPH_Caudate_FactorizedScore)

qqnorm(residuals(MPH_Caudate_FactorizedScore))

## effects of SUL vs PBO on factorized model


SUL_VS_FactorizedScore <- lmer(Score ~ 1 + Drug*VS_ki*ScoreType + Session + (1 | ID), data = df_SUL_Factorized, REML=F) 
print(summary(SUL_VS_FactorizedScore), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Score ~ 1 + Drug * VS_ki * ScoreType + Session + (1 | ID)
##    Data: df_SUL_Factorized
## 
##      AIC      BIC   logLik deviance df.resid 
##    972.1   1015.3   -475.1    950.1      365 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7721 -0.6374 -0.1463  0.4634  4.2719 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.1462   0.3824  
##  Residual             0.6207   0.7878  
## Number of obs: 376, groups:  ID, 94
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)            -0.58791    0.44416  -1.324
## Drug1                  -0.25778    0.31422  -0.820
## VS_ki                  24.68895   29.91420   0.825
## ScoreType1             -0.03210    0.31418  -0.102
## Session                 0.10353    0.05470   1.893
## Drug1:VS_ki            17.51553   21.39430   0.819
## Drug1:ScoreType1       -0.05811    0.31418  -0.185
## VS_ki:ScoreType1       33.65525   21.39408   1.573
## Drug1:VS_ki:ScoreType1  3.46534   21.39408   0.162
plot_model(SUL_VS_FactorizedScore, show.values=TRUE, value.offset = .3)

plot_model(SUL_VS_FactorizedScore,  type = "pred", terms = c("VS_ki","ScoreType","Drug"))

plot(SUL_VS_FactorizedScore)

qqnorm(residuals(SUL_VS_FactorizedScore))

SUL_Putamen_FactorizedScore <- lmer(Score ~ 1 + Drug*Putamen_ki*ScoreType + Session + (1 | ID), data = df_SUL_Factorized, REML=F) 
print(summary(SUL_Putamen_FactorizedScore), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Score ~ 1 + Drug * Putamen_ki * ScoreType + Session + (1 | ID)
##    Data: df_SUL_Factorized
## 
##      AIC      BIC   logLik deviance df.resid 
##    975.8   1019.1   -476.9    953.8      365 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8211 -0.6335 -0.1333  0.4659  4.2397 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.1463   0.3825  
##  Residual             0.6276   0.7922  
## Number of obs: 376, groups:  ID, 94
## 
## Fixed effects:
##                               Estimate Std. Error t value
## (Intercept)                 -0.1788738  0.4399031  -0.407
## Drug1                       -0.0659740  0.3121581  -0.211
## Putamen_ki                  -3.5193419 25.6718981  -0.137
## ScoreType1                   0.3914818  0.3119329   1.255
## Session                      0.1083277  0.0550036   1.969
## Drug1:Putamen_ki             3.7289515 18.4033478   0.203
## Drug1:ScoreType1            -0.0009515  0.3119329  -0.003
## Putamen_ki:ScoreType1        3.9571244 18.4008205   0.215
## Drug1:Putamen_ki:ScoreType1 -0.3987063 18.4008205  -0.022
plot_model(SUL_Putamen_FactorizedScore, show.values=TRUE, value.offset = .3)

plot_model(SUL_Putamen_FactorizedScore,  type = "pred", terms = c("Putamen_ki","ScoreType","Drug"))

plot(SUL_Putamen_FactorizedScore)

qqnorm(residuals(SUL_Putamen_FactorizedScore))

SUL_Caudate_FactorizedScore <- lmer(Score ~ 1 + Drug*Caudate_ki*ScoreType + Session + (1 | ID), data = df_SUL_Factorized, REML=F) 
print(summary(SUL_Caudate_FactorizedScore), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Score ~ 1 + Drug * Caudate_ki * ScoreType + Session + (1 | ID)
##    Data: df_SUL_Factorized
## 
##      AIC      BIC   logLik deviance df.resid 
##    975.2   1018.4   -476.6    953.2      365 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9221 -0.6350 -0.1389  0.4458  4.2110 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.1464   0.3826  
##  Residual             0.6264   0.7914  
## Number of obs: 376, groups:  ID, 94
## 
## Fixed effects:
##                             Estimate Std. Error t value
## (Intercept)                 -0.09991    0.40323  -0.248
## Drug1                       -0.15658    0.28257  -0.554
## Caudate_ki                  -9.66734   27.83538  -0.347
## ScoreType1                   0.37966    0.28246   1.344
## Session                      0.10688    0.05489   1.947
## Drug1:Caudate_ki            10.97217   20.01212   0.548
## Drug1:ScoreType1            -0.14891    0.28246  -0.527
## Caudate_ki:ScoreType1        5.60054   19.98541   0.280
## Drug1:Caudate_ki:ScoreType1 10.10048   19.98541   0.505
plot_model(SUL_Caudate_FactorizedScore, show.values=TRUE, value.offset = .3)

plot_model(SUL_Caudate_FactorizedScore,  type = "pred", terms = c("Caudate_ki","ScoreType","Drug"))

plot(SUL_Caudate_FactorizedScore)

qqnorm(residuals(SUL_Caudate_FactorizedScore))

## effects of Ki during placebo on factorized model

OnlyPBO <- dplyr::filter(FactorizedScores, Drug %in% c("PBO"))  

PBO_VS_FactorizedScore <- lmer(Score ~ 1 + VS_ki*ScoreType + Session + (1 | ID), data = OnlyPBO, REML=F) 
## boundary (singular) fit: see ?isSingular
print(summary(PBO_VS_FactorizedScore), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Score ~ 1 + VS_ki * ScoreType + Session + (1 | ID)
##    Data: OnlyPBO
## 
##      AIC      BIC   logLik deviance df.resid 
##    503.5    526.2   -244.8    489.5      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7456 -0.6122 -0.1823  0.4128  4.3174 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.0000   0.0000  
##  Residual             0.7913   0.8896  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)      -0.25688    0.51192  -0.502
## VS_ki             9.62638   34.33228   0.280
## ScoreType1        0.02601    0.50169   0.052
## Session           0.04565    0.08047   0.567
## VS_ki:ScoreType1 30.18991   34.16245   0.884
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
plot_model(PBO_VS_FactorizedScore, show.values=TRUE, value.offset = .3)

plot_model(PBO_VS_FactorizedScore,  type = "pred", terms = c("VS_ki","ScoreType"))

plot(PBO_VS_FactorizedScore)

qqnorm(residuals(PBO_VS_FactorizedScore))

PBO_Putamen_FactorizedScore <- lmer(Score ~ 1 + Putamen_ki*ScoreType + Session + (1 | ID), data = OnlyPBO, REML=F) 
## boundary (singular) fit: see ?isSingular
print(summary(PBO_Putamen_FactorizedScore), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Score ~ 1 + Putamen_ki * ScoreType + Session + (1 | ID)
##    Data: OnlyPBO
## 
##      AIC      BIC   logLik deviance df.resid 
##    504.3    527.0   -245.2    490.3      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6556 -0.6098 -0.1839  0.4319  4.2902 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.0000   0.0000  
##  Residual             0.7947   0.8915  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)           -0.04644    0.50472  -0.092
## Putamen_ki            -4.60570   29.50633  -0.156
## ScoreType1             0.39243    0.49640   0.791
## Session                0.04945    0.08086   0.612
## Putamen_ki:ScoreType1  4.35583   29.28233   0.149
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
plot_model(PBO_Putamen_FactorizedScore, show.values=TRUE, value.offset = .3)

plot_model(PBO_Putamen_FactorizedScore,  type = "pred", terms = c("Putamen_ki","ScoreType"))

plot(PBO_Putamen_FactorizedScore)

qqnorm(residuals(PBO_Putamen_FactorizedScore))

PBO_Caudate_FactorizedScore <- lmer(Score ~ 1 + Caudate_ki*ScoreType + Session + (1 | ID), data = OnlyPBO, REML=F) 
## boundary (singular) fit: see ?isSingular
print(summary(PBO_Caudate_FactorizedScore), corr=F) # print summary without fixed effect correlation matrix
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Score ~ 1 + Caudate_ki * ScoreType + Session + (1 | ID)
##    Data: OnlyPBO
## 
##      AIC      BIC   logLik deviance df.resid 
##    503.9    526.6   -245.0    489.9      181 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6763 -0.6183 -0.2027  0.3962  4.2602 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.0000   0.0000  
##  Residual             0.7931   0.8906  
## Number of obs: 188, groups:  ID, 94
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)             0.16080    0.47189   0.341
## Caudate_ki            -20.25790   31.80862  -0.637
## ScoreType1              0.52857    0.44951   1.176
## Session                 0.04874    0.08017   0.608
## Caudate_ki:ScoreType1  -4.49994   31.80427  -0.141
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
plot_model(PBO_Caudate_FactorizedScore, show.values=TRUE, value.offset = .3)

plot_model(PBO_Caudate_FactorizedScore,  type = "pred", terms = c("Caudate_ki","ScoreType"))

plot(PBO_Caudate_FactorizedScore)

qqnorm(residuals(PBO_Caudate_FactorizedScore))

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.